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Magnetic resonance imaging method 



The invention relates to a magnetic resonance (MR) method for forming an 
image of an object wherein a set of non-linear trajectories in k-space is acquired, whereas the 
density of said set of trajectories is substantially lower than the density corresponding to the 
object size. Signals along these trajectories are sampled by means of one or more receiving 
5 antennae, and a magnetic resonance image is derived from these signals and on the basis of 
the spatial sensitivity profile of the set of receiving antennae. The invention notably pertains 
to a magnetic resonance imaging method in which magnetic resonance signals are acquired 
by means of a receiver antennae system and a magnetic resonance image is reconstructed on 
the basis of the magnetic resonance signals. 
10 Such a magnetic resonance imaging method is known from the international 

application WO 01/73463. 

In this known magnetic resonance imaging method the magnetic resonance 
signals are acquired by scanning along a trajectory in k-space. The known magnetic 
resonance imaging method offers a high degree of freedom in choosing the acquisition 
15 trajectory to be followed through k-space. Notably, acquisition trajectories, notably spiral 
shaped trajectories, which give rise to irregular sampling patterns in k-space may be used. 

The invention also relates to an MR apparatus and a computer program 
product for carrying out such a method. 

Normally, in parallel imaging as SENSE (Pruessmann) or SMASH 
20 (Sodickson) the reconstruction of the image is performed by a Cartesian gridding of k-space 
or image space, respectively. 

In US -A-2003/0 122545 a magnetic resonance imaging method is described 
wherein the degree of sub-sampling is chosen such that the ensuing acquisition time for 
receiving (echo) series of magnetic resonance signals due to an individual RF excitation is 
25 shorter than the decay time of the MR signals. Preferably, a segmented scan of the k-space is 
performed, the number of segments and the number of lines scanned in each segment being 
adjustable and a predetermined total number of lines being scanned. A small number of 
segments is used such that the acquisition time for receiving the magnetic resonance signals 
for the complete magnetic resonance image is shorter than the process time of the dynamic 
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process involved. Although the method is described w^ith a scanning trajectory of straight 
lines in k-space also other trajectories like curved lines as arcs of circle or spirals could be 
possible. However, in such case more complex frequency and phase encoding of the 
magnetic resonance signals will be required. A specific solution for continuous non-Cartesian 
5 trajectories in k-space is not described. 

In an article of M. Bydder et. al. in Magn. Reson. Med 10 (2002) it is 
mentioned that partially parallel imaging techniques for reconstruction of under-sampled k- 
space data from multiple coils may be used with arbitrary acquisition schemes (e.g. 
Cartesian, spiral etc.) by casting the problem as a large linear system of equations. For 
10 realistic applications, however, the computational costs for solving this system directly is 
prohibitive. To date there is no realistic solution for a procedure of fast reconstruction of 
grossly undersampled data on a continuous non-Cartesian trajectory especially as spiral 
sampling in k-space. 

An object of the present invention is to further reduce the computational effort 
15 involved in the reconstruction of the magnetic resonance imaging method from the acquired 
magnetic resonance signals. 

■ This object is achieved by the magnetic resonance imaging method of the 
invention, wherein 

magnetic resonance signals are acquired by means of a receiver antennae system via a 
20 plurality of signal channels 

which receiver antennae system has a sensitivity profile 
the magnetic resonance signals are acquired with undersampling 

for respective orientated sector shaped regions in k-space, regularly re-sampled 
magnetic resonance signals are re-sampled on a regular sampling grid from the 
25 undersampled acquired magnetic resonance signals 

- the re-sampling involving convolution of the undersampled acquired magnetic 
resonance signals by a gridding kernel 
the gridding kernel depending on 
the orientation of the sector shaped region at issue and 
30 - the sensitivity profile of the receiver antennae system and 

a magnetic resonance image is reconstructed from the re-sampled magnetic resonance 
signals. 

The present invention is based on the following insights. In order to achieve a 
fast reconstruction, e.g. techniques such as fast Fourier transformation (FFT) are employed. 



wo 2005/078470 



PCT/IB2005/050458 



As input, these techniques require that data are sampled on a regular sampling grid. Further, a 
wide class of acquisition trajectories in k-space, notably spiral shaped trajectories and 
trajectories that include spiral segments are accurately or at least fairly approximated by 
(almost) parallel segments of the trajectory in respective sector shaped regions of k-space. In 
5 general, sector shaped regions may be regions of k-space which have a main axis that passes 
through the origin of k-space. Such sector shaped regions extend between angular 
boundaries, that is between a respective minimum and maximum modulus of the k- vector to 
the periphery of k-space and are bounded by radial boundaries that extend radially from the 
origin of k-space. The sector shaped regions maybe full sectors which extend from or through 

lO the origin of k-space into the periphery of k-space. In two dimensions the sector shaped 

regions are flat sectors or sector segments or sectors that extend point-symmetrically through 
the origin of k-space, in three dimensions the sector shaped regions are cones or portions of 
cones in k-space. According to the invention, the reconstruction which involves a re-gridding 
to re-sample the acquired magnetic resonance signals to re-sampled magnetic resonance 

15 signals on the regular grid is performed separately for the individual sector shaped regions. 
The re-gridding involves a convolution with a gridding kernel. The gridding kernel depends 
on the orientation of the sector shaped region at issue so as to accoimt for the appropriate 
direction in the image space into which aliasing will occur due to the Fourier relationship 
between pixel- values of the magnetic resonance image and the re-sampled magnetic 

20 resonance signals in k-space. Further, the gridding kernel involves the sensitivity profile of 
the receiver antermae system in order to take account of aliasing that is caused by 
undersampling of the acquired magnetic resonance signals. To derive the gridding kernel on 
the basis of the orientation of the sector shaped region and on the basis of the sensitivity 
profile does not require much computational effort. The computational effort is notably 

25 reduced because aliasing is caused by a comparatively small number of pixels or voxels. 
Accordingly, matrix inversions are only needed for matrices having a relatively low 
dimensionality. The actual re-sampling onto the regular grid, such as a Cartesian square 
lattice, involves only convolution with the gridding kernel which takes only little 
computational effort. The final reconstruction of the magnetic resonance image is then 

30 performed by a FFT technique that takes only a short computation time. 

It appears that the dependence on the orientation of the sector shaped region in 
k-space of the gridding kernel involves a quite smooth variation. Accordingly, the gridding 
kernel for a particular orientation is also accurately valid for rather wide sectors in k-space. 
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It is also an object of the present invention to provide a magnetic resonance 
imaging method enabling a fast reconstruction of grossly undersampled non-Cartesian 
sampling in k-space, especially along a spiral trajectory. It is a fixrther object of the present 
invention to provide a system and a computer program product for performing such a 
5 method. 

This object is achieved by means of a magnetic resonance imaging method 
according to the invention as claimed in particular in Claims 1, 2 and 3. 

These objects are achieved by a method as claimed in claim 1, by an MR 
apparatus as claimed in claim '6and by a computer program product as claimed in claim 21. 

10 It is a main advantage of the present invention that formulations of SPACE- 

RIP and non-Cartesian SENSE are derived that represent coil sensitivity information in the 
Fourier domain. Due to the small number of Fourier terms required, the linear system is 
highly sparse and so allows efficient solution of the equations. Thus, spiral scanning is made 
feasible at a high degree of undersampling so that a very fast acquisition and reconstruction is 

15 achieved. 

The main aspect of the present invention is that a non-Cartesian trajectory in 
k-space can be described locally by a coordinate system of imaginary parallel tangential lines 
w^hich form locally a Cartesian grid in order to perform subsampling like SENSE or SMASH. 
If the whole k-space is subdivided by rays divided homogeneously over an angle of 360° a 
20 continuous system of local Cartesian grids is obtained. These parts of k-space are than locally 
reconstructed and converted as a whole to an image. 



This and other advantages of the invention are disclosed in the dependent 
25 claims and in the following description in which an exemplified embodiment of the invention 
is described with respect to the accompanying drawings. Therein shows: 

Fig. 1 an undersampled spiral trajectory in k-space, 

Fig. 2 the same spiral trajectory as in Fig. 1 with hypothetical parallel scan 
lines for a region around the radius with an angle 0, 
30 Fig. 3 folding points in the image corresponding to a hypothetical Cartesian 

sampling pattern as shown in Fig. 2, 

Fig. 4 an apparatus for carrying out the method in accordance with the present 

invention, and 

Fig. 5 a circuit diagrgmi of the apparatus as shown in Fig. 4. 
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Basis of the present invention 

In Fig. 1 a spiral scan trajectory has been depicted, which is a single spiral arm 
2 in single shot EPI. The dots 3 represent a Cartesian grid of a density that would be required 
5 to properly image the Field-of-View (FOV) encompassing the object to be imaged. The 
actual density may also be slightly higher (so called "overgridding"), corresponding to a 
region that is slightly larger than the object The spiral arm 2 has been grossly undersampled 
according to the SENSE method with an undersampling factor of about two. This can 
immediately be seen from Fig. 1 as the distance between the spiral parts of the arm is at about 

10 a distance of two dots 3. From the point of view of the Nyquist criterion this is insufficient 
sampling. However, if the signal of that trajectory has been sampled by at least two receiver 
antennae or coils having different spatial-sensitivity patterns, the image can be reconstructed 
nevertheless. The reconstruction of the image requires the solution of a set of approximately 
N equations with N unknowns, where N is of the order of magnitude of the number of sample 

15 points times the number of coils, or of the order of the number of pixels in the resulting 
image (i.e., N is about 10"^ to 10^). This means that a number of N equations should be 
solved, which is not feasible by direct matrix inversion. Therefore, an iterative solution is 
proposed by several authors, which requires about ten iterations, each involving expensive 
computational gridding operations. 

20 

Outline of the algorithm according to the present invention 

In Figure 2 in the spiral arm 2 a radial line 5 with an angle 0 is depicted, which 
traverses the spiral arm 2. At the crossing points between the radial line 5 and the spiral arm 
2 tangential lines 6a, 6b, 6c and 6d are drawn, which show that the neighboring parts of spiral 

25 arm 2 are more or less parallel and equidistant. At present an Archimedic spiral is shown, but 
also other spiral functions may be used. This situation is well known form the Cartesian 
approach in parallel imaging: if these equidistant lines 6a to 6b would have covered the 
whole k-space, then reconstruction would be much less laborious. In image space a discrete 
number of object-points would "simply" be folded onto each other, as shown just for 

30 simplicity with two points in Figure 3. 

In principle, the problem can be solved then by "normal" SENSE 
reconstruction. This can be written as sum of receiver antenna signals irif^ (X, Y) "weighted" 

to a function a^^ {X^Y). This can also be written in the Fourier-domain as 
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piX,Y)== ^F^^{akik^,ky)®jUkik^,ky)}, or (1) 

coils k 

pi.X,Y) = F-\ '^^k(J^:„ky)®lXk(k^,ky)} (2) 

coils k 

with Mki^xiky) the measured data along the hypothetical equidistant lines 6a to 6b, with a 

the Fourier transform of (X, Y) . 

5 It is noted that equations (1) or (2) describe exactly the sanie operations 

performed for normal spiral imaging without undersampling (SENSE, SMASH). There, the 
meaning of a^{X,Y) is the "gridding kernel", which is in essence the Fourier transform of a 

box (but tapered with smooth edges to prevent that (X, Y) having a large support). 

In the present case, the shape of (X, Y) is not a tapered box, but a 

10 "reconstructing function", which depends essentially on the coil sensitivity pattern of all 
receiver antennae or coils, on the folding distance of the SENSE method and eventually 
partly on the object shape (due to regularization). Yet, since the coil sensitivity functions are 
expected to be smooth functions in space, the functions (X, Y) are also expected to be 

smooth in space. For that reason, the gridding function a^^ (k^ , A:^) is expected to have a 

15 relatively small support. It is supposed that a support of 12* 12 to 16* 16 Cartesian points 
will be sufficient (where for gridding of normal imaging a support of 4 * 4 is usually 
enough). 

The obtained gridding function a^ik^.ky) can be applied perfectly to 

reconstruct data from a set of parallel equidistant lines that are angulated with respect to the 
20 required grid. However, in this case the data is sampled along a spiral arm, and not along a 
line. That means that the obtained gridding kernel is only valid for points that are strictly 
positioned on the radius with an angle 9. Strictly the gridding kernel (X, Y) should be 

calculated for an infinity of situations. Yet, coil sensitivity patterns are normally smooth 
functions of space. This means that the weighting function a^(A^,7) (and consequently the 

25 gridding function 0Ck(kj^^ky)) will not alter significantly if the folding direction is slightly 

changed. The "folding direction" is defined by the angle between the line of the folding 
points, or, equivalently, by the orientation of the hypothetical parallel lines 6a to 6d. For that 
reason, the obtained gridding function can be applied in a predetermined region around the 
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radius with angle 0. This allows to calculate (Xkik^,ky)fox a limited number of radii (e.g. 10 
or 20). 

Resulting algorithm 

5 It is assumed that coil sensitivities are known in the entire relevant region, and 

that there is some knowledge on the presence of the object (as in Cartesian SENSE). Given is 
a spiral trajectory in two-dimensional k-space. The only relevant issue is then the distance 
between the spiral arms. Reconstruction according to the present invention will be performed 
by following steps: 

10 1. A Cartesian grid is chosen. The distance (or density) should correspond to the size of 

the object under study. 

2. A number of equidistant radii over the k-space is defined (e.g. 10 or 20). 

3. For each radius, the direction of the lines 6a to 6b tangential to the spiral arms is 
determined, whereas the distance between the tangential lines should be independent 

15 of the angle. 

4. For the obtained set of folding points, the SENSE reconstruction matrix 

^ = (S^ • • iS + )~^ • • is calculated, wherein S is the receiver antenna 
or coil sensitivity matrix, ^ is the noise covariance between the coil channels, R is the 
regularization matrix and S*' means the hermitian conjugate of S. The segments of the 
20 matrix A are combined into the function a^^ ( A^, Y) for each receiver antenna or coil k. 

(It is noted that this step is part of the normal SENSE reconstruction.) 

5. Outside of the image area, the function af^{X^Y) is set to zero (zero -padding); a 

Fourier transformation of this function is then performed into cCf^ik^^ky) , or any 

other suitable method to interpolate values of a^^ik^^ky) on a sufficiently fine grid 

25 for the subsequent convolution and resampling is used. This function is expected to be 

relevant only for small values of kx and ky, so that the outer parts thereof can be 
discarded. 

6. The data is acquired. 

7. The sampling density compensation is performed. 

30 8. For each point along the spiral trajectory, the radius that comes closest to that sample 

point is determined 



> 
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9. Using the gridding function cCkik^^ky) calculated for the closest radius of each coil, 

and the acquired data sample of each coil, a gridding operation is performed. This 
operation is part of the normal spiral scan reconstruction procedure, only the extent of 
may be larger. 

5 10. The data is summed over the coil elements. 

1 1 . The Cartesian grid points are Fourier-transformed. 

It is noted that for dynamic scans (or any type of scans in which a multitude of 
data sets for the same geometric positions is acquired), steps 1 to 5 have to done only once. 

10 Refinements of the algorithm 

Rl . A grid can be chosen which does not correspond to the normal FOV (or "the size of 

the object under study") but to a slightly larger FOV, and thus being slightly denser. 

The margins between the user selected FOV and the larger or "overgridded" FOV are 

known to contain no object. Towards the edge of the overgridded FOV, the 
15 regularization values R are gradually forced to zero. In such a manner, discontinuities 

of the functions ( , Y) are avoided, leading to a smaller support of aj^ {k^ ^^y)* 

R2. The functions (X, Y) can be preconditioned by first multiplying them with a 

common shaping function, in order to reduce the support of a^^ (^x^^y) • This may be 

a tapering window function, or a multiplication by e.g. the sum of squares of 
20 sensitivities, to prevent huge values on points in space where all coils are insensitive. 

R3. The two nearest radii can be taken and the gridding kernel functions between them 
can be interpolated. In a more efficient way, both radii are gridded and the result 
thereof is interpolated. 
R4. The most central region of k-space can be reconstructed by another method (e.g. 
25 direct inversion) and the gridding result is blended with the alternative reconstruction 

of the central region. 

R5. The functions a^^ (X, Y) are divided into a defined number of subfunctions, for which 
the support of the corresponding cci^(kx,ky)'is particularly small, and the sum of 
which equals or approximates the original af^(X,Y) in the full Field-of-View. This 
30 allows to cope with sharp transitions in (X, Y) , which most notably occur at the 

edges of the reduced Field-of- Views. A separate gridding for each set of subfunctions 
is required in this case. Compared with expanding the size of the local convolution 
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kernel adequately to represent sharp transitions, this may be computationally 
attractive. One natural choice would be to take each reduced Field-of-View 
separately, and treat the function - with or without periodic replication - with some 
sort of tapering window function. 
5 R6. For a similar reason, sets of samples assigned to adjacent radii can be separately 

gridded and transformed. This allows to provide "space" for errors with the full FOV, 
very much like the conventional oversampling provides "space" for errors at the edges 
of the full FOV. 

R7. The proposed algorithm may be used to generate an image with which one of the 
10 known iterative reconstruction algorithms (e.g. the conjugate gradient method) is then 

initialized. This allows to substantially reduce the number of iterations required by 
these methods to achieve an adequate image quality. Equally, this strategy permits to 
eliminate any artifacts potentially remaining with the proposed algorithm for 
reasonable parameter settings (i.e. limited support of (Xj^(k^^ky) and limited number 

15 of radii). 



Extension of the method 

In principle, the local neighborhood of each k-space sampling point and the 
local degree of subsampling may be considered separately. In this case, steps 1 to 5 of the 
20 method according to the present invention would be performed for sets of points with similar 
local properties, which may be arbitrarily distributed in k-space. This would allow to apply 
the proposed algorithm also to, among others, variable density spiral and conventional radial 
acquisitions. 

The apparatus shown in Fig. 4 is an MR apparatus which comprises a system 
25 of four coils 5 1 for generating a steady, uniform magnetic field whose strength is of the order 
of magnitude of from some tenths of Tesla to some Tesla. The coils 51, being concentrically 
arranged relative to the z axis, may be provided on a spherical surface 52. The patient 60 to 
be examined is arranged on a table 54 which is positioned inside these coils. In order to 
produce a magnetic field which extends in the z direction and linearly varies in this direction 
30 (which field is also referred to hereinafter as the gradient field), four coils 53 as multiple 
receiver antennae are provided on the spherical surface 52. Also present are four coils 57 
which generate a gradient field which also extends (vertically) in the x direction. A magnetic 
gradient field extending in the z direction and having a gradient in the y direction 
(perpendicularly to the plane of the drawing) is generated by four coils 55 which may be 
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identical to the coils 57 but are arranged so as to be offset 90° in space with respect thereto. 
Only two of these four coils are shown here. 

Because each of the three coil systems 53, 55, and 57 for generating the 
magnetic gradient fields is symmetrically arranged relative to the spherical surface, the field 
5 strength at the center of the sphere is determined exclusively by the steady, uniform magnetic 
field of the coil 51. Also provided is an RF coil 61 which generates an essentially uniform RF 
magnetic field which extends perpendicularly to the direction of the steady, uniform 
magnetic field (i.e. perpendicularly to the z direction). The RP coil receives an RF modulated 
current from an RF generator during each RF pulse The RF coil 61 can also be used for 

10 receiving the spin resonance signals generated in the examination zone. 

As is shown in Figure 5 the MR signals received in the MR apparatus are 
amplified by a unit 70 and transposed in the baseband. The analog signal thus obtained is 
converted into a sequence of digital values by an analog-to-digital converter 71 . The analog- 
to-digital converter 71 is controlled by a control unit 69 so that it generates digital data words 

15 only during the read-out phase. The analog-to-digital converter 71 is succeeded by a Fourier 
transformation unit 72 which performs a one-dimensional Fourier transformation over the 
sequence of sampling values obtained by digitization of an MR signal, execution being so 
fast that the Fourier transformation is terminated before the next MR signal is received. 

The raw data thus produced by Fourier transformation is written into a 

20 memory 73 whose storage capacity suffices for the storage of several sets of raw data. From 
these sets of raw data a composition unit 74 generates a composite image in the described 
manner; this composite image is stored in a memory 75 whose storage capacity suffices for 
the storage of a large number of successive composite images 80. These sets of data are 
calculated for different instants, the spacing of which is preferably small in comparison with 

25 the measurement period required for the acquisition of a set of data. A reconstruction unit 76, 
performing a composition of the successive images, produces MR images from the sets of 
data thus acquired, said MR images being stored. The MR images represent the examination 
zone at the predetermined instants. The series of the MR images thus obtained from the data 
suitably reproduces the dynamic processes in the examination zone. 

30 The units 70 to 76 are controlled by the control unit 69. As denoted by the 

down wards pointing arrows, the control unit also imposes the variation in time of the 
currents in the gradient coil systems 53, 55 and 57 as well as the central frequency, the 
bandwidth and the envelope of the RF pulses generated by the RF coil 61. The memories 73 
and 75 as well as the MR image memory (not shown) in the reconstruction unit 76 can be 
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realized by way of a single memory of adequate capacity. The Fourier transformation unit 72, 
the composition unit 74 and the reconstruction unit 76 can be realized by way of a data 
processor well-suited for running a computer program according the above mentioned 
method. 



